Methods and apparatus particularly useful for blind deconvolution

ABSTRACT

A method and apparatus for controlling an equalizer receiving the output of an unknown system in order to produce a desired response for recovering the input to the system are characterized by iteratively adjusting the equalizer such that the unknown system combined with the equalizer behaves essentially as a linear system whose (t,n) taps, for some combinations of t and n, are iteratively adjusted according to the following rule: ##EQU1## where s t ,n denotes the (t,n) tap before the iteration, s&#39; t ,n denotes the (t,n) tap after the iteration, I is a preselected integer greater then or equal to one, α i  i=1,2 . . . I are preselected scalars that may vary from iteration to iteration, and p i , q i  i=1,2, . . . I are preselected non-negative integers such that p i  +q i  ≧2.

FIELD AND BACKGROUND OF THE INVENTION

The present invention relates to methods and apparatus particularlyuseful for performing a blind deconvolution operation, e.g. forcontrolling an adjustable linear filter (sometimes called an equalizer)receiving the output of an unknown system in order to produce a desiredresponse for recovering (deconvolving) the input to the system.

The problem of blind deconvolution, as illustrated in FIG. 1 of thedrawings, arises when there is an unknown system having an output whichis observed, but whose input is unobserved. In such a system it isdesired to recover (deconvolve) the input by identifying the inverse ofthe unknown system. This is done by using an adjustable linear filter(sometimes called an equalizer) which is adjusted to recover, ordeconvolve, the input up to a constant shift (delay), and possibly aconstant phase or a sign ambiguity.

The problem of blind deconvolution is of considerable interest in manyfields of engineering and applied science. For example, in datacommunication, the unobserved input is a transmitted sequence of datasymbols (bits), and the unknown system represents the distortion in thetransmission channel (e.g., a telephone line) which causes theundesirable phenomena of inter symbol interference (ISI). The equalizer(adjustable linear filter) in this case attempts to recover thetransmitted sequence (i.e., the message) by canceling, or inverting, thedistortion caused by the channel.

A similar problem occurs in image restoration or reconstruction, whereinthe unobserved input is the original image (scene), and the unknownsystem represents the imperfections in the electronic or photographicmedium which essentially causes blurring effects. The equalizer(adjustable linear filter) in this case attempts to reconstruct theoriginal image by de-blurring, that is canceling the imperfections inthe image formation system. This problem is particularly important infields such as astronomy, remote sensing, and medical imaging, where theblurs cause a severe degradation in the quality and resolution of theobserved image. Similar deconvolution problems can be found in theanalysis of geophysical signals, in underwater acoustics, and the like.

SUMMARY OF THE INVENTION

According to the present invention there is provided a method forcontrolling an equalizer to be coupled to an unknown system having aninput and an output which equalizer receives the output of the unknownsystem in order to produce a desired response for recovering the inputto said system, characterized by iteratively adjusting the equalizersuch that the unknown system combined with the equalizer behavesessentially as a linear system whose (t, n) taps, for some combinationsof t and n, are iteratively adjusted according to the following rule:##EQU2## where s_(t),n denotes the (t, n) tap before the iteration,s'_(t),n denotes the (t, n) tap after the iteration, I is a preselectedinteger greater than or equal to one, α_(i) i=1,2 . . . I arepreselected scalars that may vary from iteration to iteration, andp_(i), q_(i) i=1,2, . . . I are preselected non-negative integers suchthat p_(i) +q_(i) ≧2.

As will be described more particularly below, the present inventioninvolves a new class of iterative deconvolution methods that convergevery fast to the desired solution in which the inverse of the unknownsystem is identified, and the input is recovered. These methods arecomputationally efficient and statistically stable. They are capable ofidentifying both the magnitude and the phase of the unknown systemtransfer function, as distinguished from the known linear predictionmethods that are incapable of identifying the phase. Another importantfeature of the proposed methods is that they do not suffer fromspurious/ambiguous solutions, and they always converge to the desiredresponse. Error analysis indicates that in many cases the method of thepresent invention is superior to linear prediction methods based on theminimum phase assumption in the sense of statistical stability.

A futher important advantage of the novel method is that it is universalin the sense that it does not impose any restrictions on the nature ofthe input, and therefore can be applied to a wide range of problems.Thus, the tap index n in the above relationships may be one-dimensional,as in the data communication problem where n represents normalizedsample time. The index n may also be multi-dimensional, as in the imagereconstruction problem.

The proposed method can be used for the purpose of system identificationwhere the equalizer can be viewed as the model for the inverse of theunknown system. The proposed method can also be used for the purpose ofspeech/image enhancement and compression, where the equalizer parameterscorrespond to the speech/image parameters, e.g. the LPC parameters inspeech.

Further features and advantages of the invention will be apparent fromthe description below.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention is herein described, by way of example only, withreference to the accompanying drawings, wherein:

FIGS. 1,2 and 6 are diagrams helpful in explaining the invention;

FIGS. 3,4,5 and 7 are flow charts illustrating four algorithms inaccordance with the present invention for iteratively adjusting theequalizer in order to identify the inverse of the unknown system and torecover (deconvolve) the unknown system input from measurements of itsinput;

FIG. 8 is a block diagram illustrating a digital, quadrature amplitudemodulation (QAM) demodulator constructed in accordance with the presentinvention, and

FIG. 9 is a flow chart illustrating the blind deconvolution method forcontrolling the demodulator of FIG. 8.

MATHEMATICAL PRELIMINARIES

It is believed a discussion of mathematical preliminaries at this pointwill be helpful in understanding the description of the presentinvention following this discussion.

A. Notations:

We denote a tensor (a matrix) by an upper-case letter (e.g. X), a columntensor (a column matrix) by a lower case bold faced letter (e.g. x), anda scalar by a lower-case non-bold letter (e.g. x).

We denote by x_(n) the n^(th) element of a set of real/complex numbers.For example, x_(n) may be the n^(th) time sample or the n^(th) componentin a sequence, or it can represent the value of a sampled image atlocation n.

Σ_(n) --sum over all possible values of n

π_(l=n) ^(m) (.)_(l) --the product (.)_(n).(.)_(n+1) . . . (.)_(m)

x^(T) --transpose of a tensor x

x⁺ --conjugate-transpose of a tensor x.

x*--complex conjugate of a scalar x.

|x|--complex magnitude (absolute value) of a scalar x.

B. Linear Systems:

The relation between the input x_(t) and the output y_(t) of a linearsystem (filter) is given by (See FIG. 2): ##EQU3## where h_(t),n is the(t,n) tap of the linear system.

The linear system whose (t,n) tap is h_(t),n convolved (combined) withthe linear system whose (t,n) tap is g_(t),n is the linear system whose(t,n) tap is Σ_(k) g_(t),k^(h) _(t-k),n-k

If the linear system is shift invariant then h_(t),n is invariant of t,in which case we denote by h_(n) =h_(t),n its n^(th) tap.

The transfer function (frequency response) of a linear shift invariantsystem whose n^(th) tap is h_(n) is given by: ##EQU4##

C. Cumulants and High-Order Spectrum

The joint cumulant of the random variables x₁,x₂, . . . x_(m) is definedby: ##EQU5## where ##EQU6## is the joint derivative of order m withrespect to the variables μ₁, . . . μ_(m) ; ln(.) is the naturallogarithm (to the base of e), j=√-1 and E{.} is the statisticalexpectation of the bracketed quantity.

A general formula relating cumulants to moments of the random variablesx₁,x₂, . . . x_(m) is given in the book: D. R. Brillinger, "Time Series,Data Analysis and Theory" Holden-Day, California, 1981.

Thus, for example: ##EQU7##

For notational convenience, we denote by: ##EQU8##

Let x_(t) be a set of random variables such that cum(x_(t);x_(t-t).sbsb.1 ; . . . x_(t-t).sbsb.p-1 ; x_(t-t).sbsb.p *; . . .x_(t-t).sbsb.p+q *) is independent of t for all combinations of t₁, . .. t_(p+q) ; Then, its (p,q)-order spectrum is defined by: ##EQU9##

D. Tensors

A tensor A is an array of real/complex valued numbers whose (n,m)element is denoted by A_(nm).

The product of a tensor A by a tensor B is a tensor C=AB whose (n,m)element is: ##EQU10##

The unit tensor is a tensor whose (n,m) element is equal to 1 whenevern=m and equal to zero whenever n≠m.

The inverse A⁻¹ of a tensor A is a tensor whose product by A is equal tothe unit tensor.

The transpose A^(T) of a tensor A is a tensor whose (n,m) element isA_(mn). The conjugate-transpose A⁺ of a tensor A is a tensor whose (n,m)element is A_(mn) *.

A column tensor b is a tensor whose (n,m) element is 0 whenever m≠0. Wedenote the (n,0) element of b by b_(n).

Finally we note that if the subscripts n and m are one-dimensional thenthe tensor A and the column tensor b are a matrix and a column-matrix(vector), respectively.

ALGORITHM IMPLEMENTATION

As described earlier, the present invention involves systems forprocessing the output of an unknown system, shown as 10 in FIG. 1, andcontrolling an adjustable linear filter (equalizer), shown as 12 in FIG.1, receiving the output in order to produce a desired response forrecovering (deconvolving) the input to the system. The inventioninvolves the development of iterative algorithms in which the unknownsystem whose (t,n) taps, for some combinations of t and n, areiteratively adjusted according to the following rule: ##EQU11## wheres_(t),n denotes the (t,n) tap before the iteration, s'_(t),n denotes the(t,n) tap after the iteration, I is a preselected integer greater thenor equal to one, α_(i) i=1,2 . . . I are preselected scalars that mayvary from iteration to iteration, and p_(i), q_(i) i=1,2, . . . I arepreselected non-negative integers such that p_(i) +q_(i) ≧2.

According to preferred embodiments of the invention described below, theequalizer is adjusted such that the unknown system combined with theequalizer behaves essentially as a linear shift invariant system whosetaps are iteratively adjusted according to the following rule: ##EQU12##where s_(n) denotes the n^(th) tap before the iteration and s'_(n)denotes the n^(th) tap after the iteration.

Preferably I=1, such that:

    s'.sub.n =αs.sub.n.sup.p (s.sub.n *).sup.q           (2)

where α is a preselected scalar that may vary from iteration toiteration, and p and q are preselected non-negative integers such thatp+q≧2.

To simplify the exposition, only the procedure in (2) has beenconsidered. To further simplify the exposition, only the case where theequalizer is a linear shift invariant system has been considered.However, the indicated generalizations may also be useful. Following aredescriptions of several algorithm implementations of the methods in (2).

If the unknown system input can be modeled as a set of statisticallyindependent random variables--which is typically the case in digitalcommunications--and the equalizer is a linear shift invariant system,then the equalizer taps are iteratively adjusted such that:

    c'=αR.sup.-1 d                                       (3)

where c' is the column tensor whose n^(th) element, c'_(n), is then^(th) tap of the equalizer after the iteration, R is the tensor whose(n,m) element is:

    R.sub.nm =cum(y.sub.t-m ; y.sub.t-n *)                     (4)

where y_(t) is the given data at the input to the equalizer, and d isthe column tensor whose n^(th) element is:

    d.sub.n =cum(z.sub.t :p;z.sub.t *:q;y.sub.t-n *            (5)

where z_(t) is the output of the equalizer given by: ##EQU13## wherec_(n) is the n^(th) tap of the equalizer before the iteration, and wherethe subscripts n and m belong to a pre-specified set.

We note that the multiplication by R⁻¹ in (3) corresponds to spectralwhitening of the observed data y_(t). A flow chart of the algorithm isillustrated in FIG. 3.

The elements of d can also be computed by: ##EQU14## Equation (5)suggests to compute joint cumulants of y_(t) and z_(t) at each iterationof the algorithm. Equation (7) suggests to first compute the cumulantsof the given data set y_(t), then to perform the iterations with respectto the taps of the equalizer until convergence is accomplished, and onlythen to generate z_(t) according to (6). A flow chart of this algorithmis illustrated in FIG. 4.

We note that the subscripts n, m, t, etc. in the above equations may beone-dimensional as in e.g. speech processing and time series analysisproblems, or multi-dimensional as in e.g. image processing problems. Ifn is one-dimensional then R becomes a matrix and d becomes a columnvector.

If the unknown system input can not be modeled as a set independentrandom variables, then it is convenient to work in the frequency domain.The resulting iterative algorithm consists of iteratively adjusting thetransfer function (frequency response) of the equalizer such that:##EQU15## where C'(w) is the transfer function of the equalizer afterthe iteration, C(w) is the transfer function of the equalizer before theiteration, S_(y) ¹,0 (w) is the (1,0)-order spectrum (power spectrum) ofy_(t), S_(y) ^(p),q (w₁, w₂, . . . ,w_(p+q)) is the (p,q)-order spectrumof y_(t), and S_(a) ¹,0 (w) and S_(a) ^(p),q (w₁ w₂, . . . ,w_(p+q)) arepreselected functions which are related to the (1,0)-order spectrum andthe (p,q)-order spectrum of the unknown system input, and where w andw₁, w₂, . . . w_(p+q-1) belong to a prespecified set of frequencies. Wenote that w may be one-dimensional or multi-dimensional depending on theapplication. A flow chart of the algorithm is illustrated in FIG. 5.

In order to ensure the convergence of all the proposed algorithms, it issuggested to select α every few iterations such that

    cum(z.sub.t ;z.sub.t *)=ηt                             (9)

where z_(t) is the output of the equalizer, and η_(t) are preselectedset of numbers, e.g. η_(t) =cum(a_(t) ;a_(t) *) if it is known a-priori.

An important feature of the proposed methods is that they do not sufferfrom spurious/ambiguous solutions, and they converge to the desiredsolution regardless of initialization. Another important feature of theproposed methods is that they are capable of solving the blinddeconvolution problem regardless of the phase characteristics of theunknown system, unlike the conventional linear prediction methods thatcan solve the problem only if the unknown system is minimum-phase.Moreover, in many cases of practical interest, the proposed methods arestatistically more stable (i.e., have smaller variance and smaller ISIthan linear prediction methods).

We may consider combining several of the proposed algorithms fordifferent pairs of p,q in order to improve the robustness andstatistical stability of the outcome results e.g. as suggested by (1).We may also consider combining linear prediction methods with theproposed methods in order to improve statistical stability and yet toobtain a method that is capable of identifying phase characteristics.

In many practically interesting situations the observed data y_(t) iscontaminated by additive noise which we denote by e_(t) as illustratedin FIG. 6. If the noise cumulants/spectra are known a-priori or can bemeasured separately, it is recommended to eliminate its effect bysubtraction, that is to compute the elements of R and d used in theiterative algorithm in (3) as follows: ##EQU16##

Similarly, to eliminate the effect of additive noise in the iterativealgorithm suggested by (8), it is recommended to replace S_(y) ^(p),q(-w₁,-w₂, . . . -w_(p+q-1), w) and S_(y) ¹,0 (w) in (8) by [S_(y) ^(p),q(-w₁,-w₂, . . . -w_(p+q-1),w)-S_(e) ^(p),q (-w₁,-w₂, . . .-w_(p+q-1),w)] and [S_(y) ¹,0 (w)-S_(e) ¹,0 (w)] respectively, whereS_(e) ^(p),q (w₁,w₂, . . . w_(p+q)) and S_(e) ¹,0 (w) are, respectively,the (p,q)-order spectrum and the (1,0)-order spectrum of e_(t).

We note that if e_(t) is Gaussian then its (p,q)-order cumulant/spectrumis identically zero whenever p+q≧2, and therefore it affects only thecomputation of covariance matrix R or the (1,0)-order spectrum S_(y) ¹,0(w).

We also note that in the blind deconvolution problem, if we want tominimize the means square restoration error between the recovered dataset at the output of the equalizer and the unobserved input of theunknown system, it is recommended to perform a final iteration of thealgorithm (3) or (8) in which only the contribution of the noise to thehigh order cumulants/spectra is subtracted, but the contribution to thesecond order statistics is not subtracted. Thus, in the case where e_(t)is Gaussian--which is a very common case--in order to minimize the meansquare restoration error we should perform a final iteration of thealgorithm without any subtraction of the noise contribution.

We note that in the implementation of the algorithm in (3) the cumulantsare typically replaced by their sample estimates based on the given dataat the input and output of the equalizer. Similarly, in the algorithmgiven by (8) the (1,0)-order spectrum and the (p,q)-order spectrum arereplaced by their sample estimates.

The iterative algorithm given by (3) can be converted into a sequentialalgorithm by performing the iterations sequentially with the data suchthat the equalizer is adjusted sequentially.

Thus, assuming for simplicity that t is one dimensional, that theequalizer is a finite impulse response (FIR) filter whose vector taps isc=(c_(L).sbsb.1 c_(L).sbsb.1₊₁ . . . c_(L).sbsb.2)^(T), and that y_(t)is zero-mean, then in case p=1, q=1 we obtain the following sequentialalgorithm in which we perform one iteration per data sample:

    c.sup.(t) =c.sup.(t-1) +β.sub.t Q.sub.t y.sub.t *z.sub.t (z.sub.t *-γ.sub.t)                                          (12)

where

    c.sup.(t) =(c.sub.L.sbsb.1.sup.(t) c.sub.L.sbsb.1.sub.+1.sup.(t) . . . c.sub.L.sbsb.2.sup.(t)).sup.T

is the vector equalizer taps after t iterations, ##EQU17## where Q₀ is apreselected matrix; c.sup.(0) is a preselected vector; L₁ and L₂ arepreselected numbers; and β_(t), γ_(t), δ_(t) are preselected sequencesof numbers.

In case p=2, q=1 we obtain the following sequential algorithm:

    c.sup.(t) =c.sup.(t-1) +β.sub.t Q.sub.t y.sub.t *z.sub.t (|z.sub.t |.sup.2 -ρ.sub.t)         (15)

where ρ_(t) is preselected. A flow chart of the algorithm is illustratedin FIG. 7.

If we eliminate the multiplying matrix Q_(t) in (15) (i.e., replacingQ_(t) by the identity matrix) then under certain choice of β_(t) andρ_(t) our algorithm reduces to Godard's algorithm for blinddeconvolution as described in the paper: D. Godard, "Self recoveringequalization and carrier tracking in two-dimensional data communicationsystems," IEEE Trans. Commun. Vol. COM-28, No. 11, pp. 1867-1875, Nov.1980. The additional multiplication by Q_(t) in our algorithm maysignificantly improve convergence rate, statistical stability, andsensitivity to additive noise.

The sequential algorithms in (12) and (15) may be particularly useful intime series analysis and modeling, in speech enhancement and compressionand the like. Similar algorithms for other combinations of p and q canalso be derived. We may also consider combining the two algorithms, andperhaps other algorithms in the class in order to improve statisticalstability and robustness of the algorithm.

DESCRIPTION OF A PREFERRED EMBODIMENT

FIG. 8 is a block diagram illustrating a known all-digital quadratureamplitude modulation (QAM) demodulator (see e.g. the book: J. G.Proakis, "Digital Communications", McGrew-Hill, 2nd edition, 1989)constructed in accordance with the present invention. It includes anequalizer 80, in the form of a linear filter which is adjustable by aprocessor 81 in accordance with the method of the present invention. Inthe example illustrated in FIG. 8, and in its respective flow chart ofFIG. 9, the integration algorithm is that of equation (3) with p=2, q=1.Thus, by controlling the equalizer 80 to recover the input signalsequence to the demodulator, the illustrated system enables operation ofthe demodulator without the need for training sequences, and evenwithout the need to know a-priori the constellation of the transmitteddata symbols.

The illustrated demodulator is a system which detects transmitted datasymbols from a measured QAM analog communication signal. The demodulatorconsists of several sub-units whose function is described below.

The sampler unit 82 converts the QAM analog communication signal into adiscrete sequence x_(t) that is sampled at the BAUD rate. Note that inthis specific application t is a one-dimensional variable correspondingto the index of the time sample.

The timing control unit 86 supplies sampling delays such that thedemodulator is synchronized with the modulator.

The baseband conversion unit 83 generates a complex-valued sequencey_(t), whose real and imaginary parts are lowpass filtered versions ofthe sampled signal multiplied by a cosine wave and a sine wave,respectively, at the modulation frequency.

The equalizer unit 80 processes the output y_(t) of the basebandconversion unit through an equalizer in order to recover the transmitteddata symbols a_(t) up to a delay and a possible phase shift.

The phase recovery unit 84 adjusts the phase of the recovered symbols atthe output of the equalizer using a decision feedback loop.

The detector unit 85 processes the output of the phase recovery unit toproduce the estimated (recovered) data symbols a_(t), obtained byminimizing the magnitude of the restoration error ε_(t) using amemoryless nearest neighbor (minimum distance) decision rule.

We shall now turn to the equalizer adjustment unit 81, which constitutesour novel contribution to the demodulator system. The proposed algorithmis a special case of (3) where p=2 and q=1. A flow chart of thealgorithm is given in FIG. 9.

The relation between the input y_(t) and the output z_(t) of theequalizer is given by: ##EQU18## where the parameters of the equalizerL₁, L₂ and c_(L).sbsb.1, . . . c_(L).sbsb.2 are adjusted by theequalizer adjustment unit. Detailed description of the proposedalgorithm is given below:

(A) Algorithm Initialization

Accumulate N samples y₁, . . . y_(N) at the input to the equalizer.

Set L₁ =-M₀, L₂ =M₀ where M₀ is preselected positive integer (M₀ is onthe order of the length of the equalizer needed to recover thetransmitted data symbols).

Set the L×1 vector of equalizer taps c=(c_(L).sbsb.1, . . .c_(L).sbsb.2)^(T) such that c₀ =1 and c_(n) =0 for all n≠0, where L=L₂-L₁ +1 is the number of equalizer taps.

(B) Cumulants estimation

Compute the L×L matrix R whose (n,m) element is: ##EQU19## Compute theL×1 vector d=(d_(L).sbsb.1, . . . d_(L).sbsb.2)^(T) whose n^(th) elementis ##EQU20## Comment: If the cumulants of the additive noise at theinput to the equalizer are known a-priori can be measured separately, itis recommended to subtract them from R and d as suggested in (10) and(11) respectively.

(C) Equalizer Taps Adjustment

Compute the L×1 vector c'=(c'_(L).sbsb.1, . . . c'_(L).sbsb.2)^(T) :c'=R⁻¹ d

Where we note that this computation can be carried out efficiently(order of L² operations) using the Toepliz form of R (See the book: I.C. Gohberg and I. A. Feldman, "Convolution Equations and ProjectionMethods for their Solution", American Mathematical Society, 1974).

(D) Normalization

Compute the L×1 vector c"=(c"_(L).sbsb.1, . . . c"_(L).sbsb.2)^(T) :##EQU21## where η>0 is preselected, e.g. η=E{|a_(t) |² } is the averagepower of the input symbol if it is known a-priori.

(E) Equalizer's Domain Adjustment

Calculate ##EQU22## where μ≧0 is preselected, e.g. μ=0.01. Comment: Dmeasures the distance between the equalizer taps before and after theiteration.

Set L'₁ to be the maximal integer K such that ##EQU23## Set L'₂ to bethe minimal integer K such that ##EQU24## where 0<γ<1 is preselected(e.g. γ=0.1) and M₁, M₂ are preselected positive integers.

Set ##EQU25## Set

    L.sub.1 =L'.sub.1, L.sub.2 =L'.sub.2

(F) Filtering

Compute ##EQU26##

(G) Stop?

If D>0 then

Goto step (B)

else

Stop equalizer adjustment (exit).

After the iterations of the proposed algorithm have stopped, we mayswitch to a decision directed equalization (DDE) provided that theconstellation of the input data symbols is known a-priori or can beidentified from the recovered symbols at the output of the equalizer. Wethen switch back from the DDE mode of operation to the proposedalgorithm whenever is needed.

While the invention has been described with respect to severalalgorithms and one application, it will be appreciated that many othervariations, modifications and applications of the invention may be made.

What is claimed is:
 1. A method for controlling an equalizer to becoupled to an unknown system having an input and an output whichequalizer receives the output of the unknown system in order to producea response recovering the input to said system, characterized byiteratively adjusting the equalizer such that the unknown systemcombined with the equalizer behaves essentially as a linear system whose(t,n) taps, for some combinations of t and n, are iteratively adjustedaccording to the following rule: ##EQU27## where s_(t),n denotes the(t,n) tap before the iteration, s'_(t),n denotes the (t,n) tap after theiteration, I is a preselected integer greater than or equal to one,α_(i) i=1,2. . . I are preselected scalars that may vary from iterationto iteration, and p_(i), q_(i) i=1,2, . . . I are preselectednon-negative integers such that p_(i) +q_(i) ≧2.
 2. The method accordingto claim 1 wherein the unknown system combined with the equalizerbehaves essentially as a linear shift invariant system whose taps areinteratively adjusted according to the following rule: ##EQU28## wheres_(n) denotes the n^(th) tap before the iteration and s'_(n) denotes then^(th) tap after the iteration.
 3. The method according to claim 2,wherein I=1, such that

    s'.sub.n =αs.sub.n.sup.p (s.sub.n *).sup.q

where α is a scalar that may vary from iteration to iteration, and p andq are preselected non-negative integers such that p+q≧2.
 4. The methodaccording to claim 3, wherein the equalizer is a linear shift invariantsystem whose taps are iteratively adjusted such that:

    c'=αR.sup.-1 d

where c' is the column tensor whose n^(th) element, c'_(n) is the n^(th)tap of the equalizer after the iteration, R is the tensor whose (n,m)element is:

    R.sub.nm =cum(y.sub.t-mi ;y.sub.t-n *)

where y_(t) is the input to the equalizer, and d is the column tensorwhose n^(th) element is:

    d.sub.n =cum(z.sub.t :p; z.sub.t *q;y.sub.t-n *)

where z_(t) is the output of the equalizer given by: ##EQU29## wherec_(n) is the n^(th) tap of the equalizer before the iteration, and thesubscripts n and m belong to a pre-specified set.
 5. The methodaccording to claim 4 wherein the elements of d are computed as follows:##EQU30##
 6. The method according to claim 3 wherein the equalizer is alinear shift invariant system whose transfer function is iterativelyadjusted such that: ##EQU31## where C'(w) is the transfer function ofthe equalizer after the iteration, C(w) is the transfer function of theequalizer before the iteration, S_(y) ¹,0 (w) is the (1,0)-orderspectrum (power spectrum) of y_(t), S_(y) ^(p),q (w₁, w₂, . . . ,w_(p+q)) is the (p,q)-order spectrum of y_(t) where p and q arepreselected non-negative integers such that p+q≧2, S_(a) ¹,0 (w) andS_(a) ^(p),q (w₁, w₂, . . . , w_(p+q)) are preselected functions, andwhere w and w₁, w₂, . . . w_(p+q-1) belong to a prespecified set offrequencies.
 7. The method according to claim 3 wherein α is selected atsome iterations such that

    cum(z.sub.t ;z.sub.t *)=η.sub.t

where z_(t) is the output of the equalizer and η_(t) is preselected. 8.The method according to claim 1 wherein the cumulants/spectra ofadditive noises at some iterations are subtracted from the respectivedata cumulants/spectra.
 9. The method according to claim 1 wherein thecumulants/spectra are replaced by their estimates based on the availabledata at the input and output of the equalizer.
 10. The method accordingto claim 1 wherein the iterations are performed sequentially such thatthe equalizer is adjusted sequentially.
 11. The method according toclaim 10 wherein t is one-dimensional, and the equalizer is adjustedsequentially as follows: ##EQU32## is the vector equalizer taps after titerations, ##EQU33## where Q₀ is a preselected matrix; c.sup.(0) is apreselected vector; L₁ and L₂ are preselected numbers; and β_(t), γ_(t),δ_(t) are preselected sequences of numbers.
 12. The method according toclaim 11 wherein

    c.sup.(t) =c.sup.(t-1) +β.sub.t Q.sub.t y.sub.t *z.sub.t (|z.sub.t |.sup.2 -ρ.sub.t)

where ρ_(t) is preselected.
 13. Apparatus to be coupled to an unknownsystem having an input and an output for producing a response recoveringthe input to said system, comprising an equalizer receiving the outputof the unknown system and a processor for iteratively adjusting theequalizer such that the unknown system combined with the equalizerbehaves essentially as a linear system whose (t, n) taps, for somecombinations of t and n, are iteratively adjusted according to thefollowing rule: ##EQU34## where s_(t),n denotes the (t, n) tap beforethe iteration, s'_(t),n denotes the (t, n) tap after the iteration, I isa preselected integer greater than or equal to one, α_(i) i=1,2 . . . Iare preselected scalars that may vary from iteration to iteration, andp_(i), q_(i) i=1,2, . . . I are preselected non-negative integers suchthat p_(i) +q_(i) ≧2.
 14. Apparatus according to Claim 13 wherein theunknown system combined with the equalizer behaves essentially as alinear shift invariant system whose taps are iteratively adjustedaccording to the following rule: ##EQU35## where s_(n) denotes then^(th) tap before the iteration and s'_(n) denotes the n^(th) tap afterthe iteration.
 15. Apparatus according to claim 14, wherein I=1, suchthat

    s'.sub.n =αs.sub.n.sup.p (s.sub.n *).sup.q

where α is a scalar that may vary from iteration to iteration, and p andq are preselected non-negative integers such that p+q≧2.
 16. Apparatusaccording to claim 15, wherein the equalizer is a linear shift invariantsystem whose taps are iteratively adjusted such that:

    c'=αR.sup.-1 d

where c' is the column tensor whose n^(th) element, c'_(n) is the n^(th)tap of the equalizer after the iteration, R is the tensor whose (n,m)element is:

    R.sub.nm =cum(y.sub.t-m ;y.sub.t-n *)

where y_(t) is the input to the equalizer, and d is the column tensorwhose n^(th) element is:

    d.sub.n =cum(z.sub.t :p;z.sub.t *:q;y.sub.t-n *)

where z_(t) is the output of the equalizer given by: ##EQU36## wherec_(n) is the n^(th) tap of the equalizer before the iteration, and thesubscripts n and m belong to a pre-specified set.
 17. Apparatusaccording to claim 16 wherein the elements of d are computed as follows:##EQU37##
 18. Apparatus according to claim 15 wherein the equalizer is alinear shift invariant system whose transfer function is iterativelyadjusted such that: ##EQU38## where C'(w) is the transfer function ofthe equalizer after the iteration, C(w) is the transfer function of theequalizer before the iteration, S_(y) ¹,0 (w) is the (1,0)-orderspectrum (power spectrum) of y_(t), S_(y) ^(p),q (w₁, w₂, . . . ,w_(p+q)) is the (p,q)-order spectrum of y_(t) where p and q arepreselected non-negative integers such that p+q≧2, S_(a) ¹,0 (w) andS_(a) ^(p),q (w₁, w₂, . . . , w_(p+q)) are preselected functions, andwhere w and w₁, w₂, . . . w_(p+q-1) belong to a prespecified set offrequencies.
 19. Apparatus according to claim 15 wherein α is selectedat some iterations such that

    cum(z.sub.t ;z.sub.t *)=η.sub.t

where z_(t) is the output of the equalizer and η_(t) is preselected.